module adv_case1_mod

  use const_mod

  implicit none

  private

  public adv_case1

  real(8), parameter :: radius  = 6371.220d3
  real(8), parameter :: T0      = 12.0d0 * 86400.0d0
  real(8), parameter :: u0      = 2.0d0 * pi * radius / T0
  real(8), parameter :: R0      = radius / 3.0d0
  real(8), parameter :: alpha   = 0.5d0 * pi
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: period  = 12 * 86400

contains

  subroutine adv_case1(time_in_seconds, lon, lat, u, v, gh, ghs, g, radi)
    real(8), intent(in ) :: time_in_seconds
    real(8), intent(in ) :: lon
    real(8), intent(in ) :: lat
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: gh
    real(8), intent(out), optional :: ghs
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: radi

    real(8) lon1, lat1, lon2, lat2, r1, r2, r, b, c
    real(8) pi2

    real(8) k, cos_T
    real(8) c1, c2
    
    if(present(radi))radi = radius
    
    pi2  = pi * 2
    lon1 = pi * 5.0d0 / 6.0d0
    lat1 = 0.0d0
    lon2 = pi * 7.0d0 / 6.0d0
    lat2 = 0.0d0
    
    r = 0.5d0 * radius
    b = 0.1d0 * gravity
    c = 1.0d0 * gravity
    
    r1 = calc_distance(lon1, lat1, lon, lat)
    r2 = calc_distance(lon2, lat2, lon, lat)
    if ((r1 <= r .and. abs(lon - lon1) >= r / 6.0d0 / radius) .or. (r2 <= r .and. abs(lon - lon2) >= r / 6.0d0 / radius)) then
      if(present(gh)) gh = c
    else if (r1 <= r .and. abs(lon - lon1) < r / 6.0d0 / radius .and. lat - lat1 < -5.0d0 / 12.0d0 * r / radius) then
      if(present(gh)) gh = c
    else if (r2 <= r .and. abs(lon - lon2) < r / 6.0d0 / radius .and. lat - lat2 >  5.0d0 / 12.0d0 * r / radius) then
      if(present(gh)) gh = c
    else
      if(present(gh)) gh = b
    end if
    
    k = 10.0d0 * radius / period
    c1 = pi2 * time_in_seconds / period
    c2 = pi2 * radius / period
    cos_t = cos(pi * time_in_seconds / period)
    
    if (present(u  )) u   = k * sin(lon-c1)**2  * sin(2 * lat) * cos_t + c2 * cos(lat)
    if (present(v  )) v   = k * sin(2 * (lon-c1)) * cos(lat) * cos_t
    if (present(g  )) g   = gravity
    if (present(ghs)) ghs = 0

  end subroutine adv_case1

  ! spherical distance on unit sphere

  real(8) function calc_distance(lon1, lat1, lon2, lat2) result(res)

    real(8), intent(in) :: lon1
    real(8), intent(in) :: lat1
    real(8), intent(in) :: lon2
    real(8), intent(in) :: lat2

    res = radius * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))

  end function calc_distance
end module adv_case1_mod
